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ABSTRACT 

We investigate the gravitational collapse of rapidly rotating relativistic supermassive stars by means 
of a 3+1 hydrodynamical simulations in conformally flat spacetime of general relativity. We study the 
evolution of differentially rotating supermassive stars of q = J/M 2 ~ 1 ( J is the angular momentum 
and M is the gravitational mass of the star) from the onset of radial instability at R/M ~ 65 (R is 
the circumferential radius of the star) to the point where the conformally flat approximation breaks 
down. We find that the collapse of the star of q > 1, a radially unstable differentially rotating star 
form a black hole of q < 1. The main reason to prevent formation of a black hole of q > 1 is that 
quite a large amount of angular momentum stays at the surface. We also find that most of the mass 
density collapses coherently to form a supermassive black hole with no appreciable disk nor bar. In 
the absence of nonaxisymmetric deformation, the collapse of differentially rotating supermassive stars 
from the onset of radial instability are the promising sources of burst and quasinormal ringing waves 
in the Laser Interferometer Space Antenna. 

Subject headings: black hole physics — gravitation — gravitational waves — hydrodynamics — insta- 
bilities — relativity — stars: rotation 
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1. INTRODUCTION 

There is increasing evidence that supermassive black 
holes (SMBHs) exist at the center of all galaxies, and 
that they are the s ources whic h power active galactic 
nuclei and quasars <)Reesl 11998ft . For example, VLBI 
observations of the Keplerian disk around an object in 
NGC4258 indicate that the central object has a mass 
M ~ 2.0 x 10 7 M Q . Also, large numbers of observations 
are provided by the Hubble space telescope suggesting 
that SMBHs exist in galaxies such as M31 (7.0 x 1O 7 M ), 
M87 (3.4 x 10 9 Mp) and our own galaxy (3.7 x 10 6 M Q ) 
(see for example. lKormendvll2004l for a brief overview). 
Although evidence of the existence of SMBHs is com- 
pelling, the act ual formati on process of these objects is 
still uncertain ijReesI l200lft . Several different scenarios 
have been proposed, some based on stellar dynamics, oth- 
ers on gas hydrodynamics, and still others which combine 
the processes. At present, there is no definitive observa- 
tion as yet which confirms or rules out any one of these 
scenarios. 

Here we discuss the collapse of a supermassive star 
(SMS) as one scenario of formation of an SMBH. This 
subject is also interesting from a viewpoint of general rel- 
ativity, since the onset of radial instability at the mass 
shedding limit takes place at J/M 2 sa 1 (J is the an- 
gular momentum, M is the tota l gravitational energy) 
in the uniformly ro tating SMS l|Baumgarte fc Shapiro! 
119991 IShibatall2004h . which is also considered as a crit- 
ical threshold between formation of a stationary black 
hole (BH) and of a naked singularity. The BH unique- 
ness theorem says that a complete gravitational collapse 
of a body always results in a BH rather than a naked 
singularity, and that the final stat e of the B H should 
go into a stationary Kerr BH (e.g., Waicfl Il984h . How- 
ever there is a candidate to except this cosmic censor. 
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The collapse of a prolate spheroid with large semima- 
jor axis leads t o spindle singularities without an ap- 
parent horizon ||N^ainura.Sha piro fc Teukolskv| 119881 
IShaniro TeukolskvUl 991L 11 992D 1 . Therefore, a collapse 
of a star with the critical value J/M 2 ~ 1 may show us 
an interesting phenomenon in general relativity. What 
happens when the star of J/M 2 ~ 1 collapses? What 
is the final fate of a collapsing star of J/M 2 ~ 1? Nu- 
merical simulations can clearly show the answers to these 
questions. 

The gravitational collapse of a star of J/M 2 ~ 1 has 
been investigated most of all in axisymmetric spacetime. 
The p ioneering study in this field was made bv lNakamural 
(1981). He set up a differentially rotating star, adding 
radial velocity to induce the collapse. He found that 
the criterion of formation of a BH apparent horizon is 
J/M 2 ss 0.95. The follow ing study has been done by 
iStark & Piranl l|1985l 0j86). They set up a uniformly 
rotating n = 1 polytropic star and deplete 99 % of 
the whole pressure to induce the collapse. They found 
that the criterion of BH formation is a cr n/M ± 0.2, and 
when the star exceeds the above value, flattened disk 
formed (a ait /M = 1.2 for the case of 99 % pressure 
deplete). Note that a cr it(= J/M) represents the criti- 
cal Ker r parameter th at the collapse proceeds to a BH. 
Finally Shibata (2000) performed a collapse of a differ- 
entially rotating n — 1 polytropic star. He also depleted 
the pressure of the star to induce the collapse. He found 
that when J/M 2 is less than 0.5, a BH is formed when 
the rest mass is larger than the maximum mass of the 
pressure depleted J constant sequence. When J/M 2 is 
slightly less than 1, a BH formed when the rest mass is 
sufficiently larger than the maximum allowed mass on the 
pressure depleted J constant sequence. He also found, 

1 In order to prove violation of cosmic censor precisely, one might 
compute a global event horizon instead of an apparent horizon, 
which is difficult to treat numerically in dealing with a singularity. 
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by comparing his results from polytropic evolution and 
from T-law evolution of the hydrodynamics, that shock 
heating prevents the prompt collapse to a BH under the 
condition of J/M 2 ~ 1. 

At the same time we have submitted our paper, several 
groups have also investigated the collapse of a relativistic 
star of J/M 2 - 1 in full general relativity. iShibatal (|2004T) 
investigated the onset of radial instability of uniformly 
rotating stars in 2D and find that the criterion of J/M 2 is 
> 1 in a very soft polytropic equat ion of state (n sa 3.01 
- 3.05 ; n is a polytropic index). iDuez. Shapiro, fc Yd 
( 2004) investigated the collapse of a differentially rotat- 
ing n = 1 polytropic star in 3D by depleting the pressure 
and find that the criterion of BH formation is J/M 2 w 1, 
and when the star exceeds the above value, the collapsing 
star for ms a torus which fragments in to nonaxisymmetric 
clumps. Sckiguchi & Sh ibatal 1)2004') studied the collaps- 
ing star of J/M 2 ~ 1 in 2D with a polytropic index of 
n = 1 and 2 and find that the criterion for BH formation 
is determined by the central q c (= j/m 2 ; j is the specific 
angular momentum and m is the cylindrical mass) at the 
initial state of the star. 

The purpose of this paper is the following threefold. 
The first is to verify the nature of the gravitational col- 
lapse of a differentially rotating radially unstable star 
from the viewpoint of cosmic censor. If we collapse a 
star of J/M 2 > 1, we expect that the star cannot di- 
rectly form a BH of J/M 2 > 1 because of cosmic censor. 
Therefore it is important to find the main cause to pre- 
vent BH formation of J/M 2 > 1. From the prev ious 
computation al results, s hock heating ( Shibata 200 0]) and 
core bounce jNakamuralll98lUStark fc Piranlll985lll986t 
IShiba ta 2000) are the dominant phenomena to prevent 
a star from forming a BH. However, all of the previous 
calculations have set up violent initial data sets, adding 
radial velocity or depleting pressure, to induce the col- 
lapse, it may cause abrupt transport of the energy. We 
therefore set up a mild, natural situation, that is a col- 
lapse of differentially rotating stars from the onset of 
radial instability, to focus on a graduate transport of the 
energy. Also, we compute the collapsing star in 3D to 
allow shock propagation or bar formation, if it occurs. 

The second is to determine the final outcome of the col- 
lapse of differentially rotating SMSs. Two different types 
of rotation profile arises during the quasi-static evolution 
of the rotating star, that depends on the environment. 
One is a uniformly rotating star, maintained uniform ro- 
tation during the quasi-static evolution by sufficiently 
strong viscosity or strong magnetic field. The other is a 
differentially rotating star, impossible to retain uniform 
rotation in low viscosity and low magnetic field (even the 
star initially takes uniform rotation) . For the collapse of 
a unif ormly rotat i ng SM S from the onset of radial insta- 
bility, [Saij^OD <|2002j) studied 3D relativistic hydrody- 
namic simulation in the post-Newtonian gravitation and 
found that the collapse is coherent and that it is likely 
to form an SMBH with no significant bar nor disk for- 
mation. Followup computation has been performed in 
2D hydrodynamics in full general relativity and found 
that approximately 10 % of the total mass can form a 
disk, while approximately 90 % of that should form a BH 
l)Shibata fc Shapiro! I2002|) . For the differentially rotat- 
ing stars, the final outcome may strongly depends on the 
amount and the distribution of initial angular momen- 



tum, and the nature of angular momentum distribution 
of SMSs is largely unknown. We also do not know what 
path does a differentially rotating SMS take during the 
quasi-static evolution. When the amount of angular mo- 
mentum of the star is sufficient, the star seems to follow 
in a quasi-static evolu tion up to the mass-shedding limit. 
iNew fc Shapiro! l)2001[) find in the Newtonian quasi-static 
evolution that bar formation is inevitable before reach- 
ing mass-shedding limit. When the amount of angular 
momentum of the star is not sufficient, the star seems to 
collapse at the onset of radial instability before reaching 
secular instability (T/W « 0.14 f or a uniformly rotatin g 
incompressible Newtonian stars (|Chandrasekharl fl969')') 
or mass-shedding limit. The purpose of this paper is to 
examine the collapse of the SMS from the onset of radial 
instability. When we consider the collapse of radially un- 
stable differentially rotating stars, the final outcome has 
a possibility of differing from the collapse of uniform rota- 
tion because of the strong centrifugal force at the central 
core, which may prevent the prompt collapse. What is 
the final fate of the collapse? Does the star fragment due 
to the growth of the degree of differential rotation? Docs 
the disk form during the collapse? Relativistic simula- 
tion can clearly answer to these questions. 

Finally, it is important to probe whether the collapse of 
differentially rotating SMSs could be a promising source 
of gravitational waves. Direct detection of gravitational 
waves by ground based and space based interferometers is 
of great importance in general relativ ity, in astrophysics, 
and in cosmology (e.g.. lThornei ri998'). The catastrophic 
collapse is o ne of the p romising sources of gravitational 
waves (e.g.,|Ne3|2Q03). For a gravitational collapse of 
the star in a dynamical timescale, there are two main 
reasons that prevent a prompt collapse of the star, which 
should produce gravitational waves at that moment. One 
reason is core bounce and/or shock heating. Suppose 
gravitational force is balanced to the centrifugal force 
(M/R 2 ~ Rti 2 ) in Newtonian gravity. The total mass 
and the angular momentum (J ~ MR 2 Sl) are almost 
conserved during the collapse, assuming that gravita- 
tional radiation takes a little role during the collapse 
in the absence of nonaxisymmetric deformation prior to 
form a BH. Note that is angular velocity, R is the 
radius. We may estimate the bounce radius of the star 
in a global sense as -Rbounce ~ M(J/M 2 ) 2 . Since the 
horizon radius is roughly the order of M , core bounce 
might take place for the case J/M 2 > 1, whic h pre- 
vents violation of cosmic censor i|Nakamuralli98l . The 
other is bar formation. From the dimensional analysis, 
we can describe the rotational kinetic energy T and the 
gravitational binding energy W as T ~ MR 2 ft 2 and 
W ~ M 2 /R. We also accept the assumption that to- 
tal mass and angular momentum are almost conserved 
during the collapse. We could estimate the radius of 
bar formation in terms of the ratio of the rotational 
kinetic energy to the gravitational binding energy as 
i?bar « (M/R)(T/Wy 1 (J/M 2 ) 2 . Since the dynamical 
instability for a uniformly rotating , incompressible New- 
tonian star sets in at T/W ~ 0.27 l|Cnandrasckhar fl969|) 
and for relativistic gravitation as T/W ~ 0.24 — 0.2 6 
( Shib ata. Baumgarte! fc Shapircl2000tlSaiio et al.l 200D , 
bar formation takes place at the radius R ~ AM for the 
collapse of J/M 2 ~ 1. Therefore, we may expect that a 
gravitational collapse of J/M 2 ~ 1 is a promising source 
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of quasi-periodic gravitational waves. 

This paper is organized as follows. In §|21we present 
basic equations of our conformally flat approximation in 
general relativity. We demonstrate our code tests in § 02 
We discuss our initial data sets and our numerical results 
for our catastrophic collapse in §0J In §Elwe summarize 
our findings. Throughout this paper, we use geometrized 
units (G = c = 1) and adopt Cartesian coordinates 
(x,y,z) with the coordinate time t. Greek and Latin 
indices run over (t, x, y, z) and (x, y, z), respectively. 

2. 3+1 RELATIVISTIC HYDRODYNAMICS IN 
CONFORMALLY FLAT SPACETIME 

In this section, we briefly describe the co nformally 
flat spacetime (e.g., Ilsenberg fc Nesteil I1980D . We 
solve the fully relativistic equations of hydrodynam- 
ics, but neglect nondiagonal sp atial metric components 
(|Wilson & Mathewslli9Mll99^ . 



2.1. The gravitational field equations 

We define the spatial projection tensor = g 
where g^ v is the spacetime metric, 



m = 



(1/a, —p) l /a) the unit normal to a spatial hypersurface, 
and where a and j3 l are the lapse and shift. Within a 
first post-Newtonian approximation, the spatial metric 
9ij — lij may always be chosen to be conformally flat 



jii = ip Sij, 



(1) 



where ip is the conformal factor (see 


Chandrasekhar 1965; 


Blanchet. Damour. & Schafcr 1989 


). The spacetime line 



element then reduces to 

ds 2 = (-a 2 + (3 k (3 k )dt 2 + 2f3 t dx l dt + ip A 5 l3 dx l dx ] . (2) 

We adopt maximal slicing, for which the trace of the 
extrinsic curvature K^j vanishes, 



K = j^Kij = 0. 



(3) 



The gravitational field equations in conformally flat 
spacetime for the five unknowns a, (3 1 , and ip can then 
be derived con veniently from the 3+1 formalism (e.g., 
ISaiio et al.ll2fiml) . 

Since the spatial metric is conformally flat, the trans- 
verse part of its time derivative vanishes. The transverse 
part of the evolution equation of the spatial metric there- 
fore relates the extrinsic curvature to the shift vector, 



2aiP~ i K lJ = 5M 1 + Sudj/3 1 - -5M 1 



(4) 



Inserting eq. (0J into the momentum constraint equa- 
tion, then, yields an equation for the shift j3 l 



0* In 



x (d^ + SuP k d k p l - Uidtf 1 



(5) 



where A = S^didj is the flat space Laplacian and Ji = 
—n^^jT^u is the momentum density. In the definition 
of Ji, T^ v is the stress energy tensor. 

The conformal factor ip is determined from the Hamil- 
tonian constraint 



Alp 



1 



(6) 



where pn = n^n^T^ is the mass-energy density mea- 
sured by a normal observer. 

Maximal slicing implies dfK = 0, so that the trace of 
the evolution equation for the intrinsic curvature yields 
an equation for the lapse a 

A(a<ip) = 2vra7/> 5 (pH + 2S*) + |cm/> 5 KijK ij , (7) 

8 

combined with eq. ©, where S = jjkT^ . 

Therefore, conformally flat gravitational field equa- 
tions for the five unknowns ip, aip, (3 % can be derived 
by eqs. © - Q. 

2.2. The matter equations 

For a perfect fluid, the energy momentum tensor takes 
the form 



T^ = p(l 



P 

E+ — 
P 



u»u v + Pg>- 



(8) 



where p is the rest-mass density, e the specific internal 
energy, P the pressure, and u 11 the four-velocity. 
We adopt a T-law equation of state in the form 

P=(r-l)pe, (9) 

where T is the adiabatic index which we set w 1.33 ~ 1.34 
in this paper. 

In the absence of thermal dissipation, eq. ©, together 
with the first law of thermodynamics, implies a poly- 
tropic equation of state 

P = K P 1 + 1/n , (10) 

where n = 1/(T — 1) is the polytropic index and k is a 
constant. 

From V p T' Ily = together with the equation of state 
(eq. 0), we can derive the energy and Euler equations 
according to 

de* d(e*v j ) I ,_i+i/r p d f t ,6 ,i Vv n 



dxi 

d(p^Ui) d(p„UiV J ) 



dt 



dx? 



-aiP 6 (P + P vis ), 



p*au a t i 



2p*u k u k 



ip 5 u 



(12) 



(13) 
(14) 
(15) 



where 

e* = (pe) 1 / r cu t V 6 , 

i dx l u l 
V = ^it" u*' 
p* = pau l ip Q , 

u t = (l+Te)u t , (16) 
u i = (l+re)u i! (17) 

and v l , P v is is the 3- velocity, pressure viscosity, respec- 
tively. Note that we treat the matter fully relativistically; 
the conformally flat approximation only enters through 
simplifications in the coupling to the gravitational fields. 
Note also that we include an artificial viscosity, since core 
bounce might occur in our gravitational collapse and that 
should produce shocks. We will explain our form of an ar- 
tificial viscosity in § O As a consequence to treat shocks 
we also need to solve the continuity equation 

dp* d(p*v l ) 
~dt + dx* 



0. 



(18) 
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separately. 

We solve the matter evolution in second order accurate 
in spa c e and time with the transfe r scheme of Ivan Leerl 
l)1977jl : IOo"hara fc Nakamural l)1989fl . 

2.3. Numerical techniques for solving gravitational field 

equations 

There are three key issues to solve gravitational field 
equations numerically. The first is to introduce symmet- 
ric tensor An related to th e extrinsic curvature as 
llShiha,ta fc Naka,murallT99! 

A 



^ 2 



A ij EE 



K ij 



->y lJ K 



(19) 



(20) 



Therefore, we describe the equations of conformal factor 
(eq. UJ) and lapse (eq. [7j) as 



7 



(21) 



A(aip) = 2iraip(pn + 25) + -aip- 7 A ij A ij . (22) 

8 

Note that the dependence value (such as ip, aip) of the 
source term is always lower than O^ 1 ), 0((a?/') 1 ), which 
is safe for a convergence during the iteration. 2 

The second is to decompose the vector type elliptic 
equation. The equation for the shift, 



1 



16naJi 



8 ij A<3 j + -didt/3 1 = I d 3 a - 6-5^J A 

= Ji, (23) 

can be further simplifie d by introduci ng a vector Pi and 
a scalar r\ according to ( Shibatal ll997l) 

AP i = J i! (24) 
Ar) = -j i x i . (25) 
The shift can then be computed from 



Sij/P = lp i -l(d i v + x k d i P k ), 



(26) 



and will automatically satisfy eq. (|23[l . Next, we decom- 
pose the symmetric tensor A y - as (e.g.. IYork1ll979|) 



At 



A* 



where 



djA*' 3 



7ijA* 



(iw)ij = diWj + djWi - -8 l3 d k w k . 



(27) 

(28) 
(29) 



Therefore, we can rewrite the Momentum constraint as 



AWi + \didjW' = 87r0 6 Ji. 



(30) 



We also decompose the vector W % with the same manner 
as the shift; 



7 



1 



SijW 3 = -Bi - -(dix + xfdiB k ) 



(31) 



2 Since we compu te th e variabilities of the matter as p*,e*,Ui, 
the first term of eq. 1211 behaves as — 2inp~ 1 (ip 6 pn). 



where Bi and \ satisfies 

AP, = 8W 6 ^, 
A X = -&mP 6 J i x t 



(32) 
(33) 



Note that Ji only appears in the presence of matter, 
which means that the source terms of eqs. (|32|l and l|33|l 
are compact. Therefore, we can compute these values 
quite accurately. 

To summarize, we have reduced Einstein equations in 
a conformally flat spacetime to 10 elliptic equations for 
10 variables (B i: x> a% Pi P; v)i 



Ax = -S^J 6 JiX 1 = 4tt5 
1 



Aip = -2Trip a pn 



A{aip) = 27ra-0( / o H + 25) + 
ee47t5 q ^, 
AP i =4naJ l ee 4ttS Pi , 
An = —AiraJiX 1 = 47r5 ?; 



X' 

7 



4tt5 



-aip 



Aij A 3 



(34) 
(35) 

(36) 

(37) 
(38) 
(39) 



These Poisson-type equations are solved by imposing the 
following boundary condition at outer boundaries 



B x 

B v 
B z 

X 

aip 
P, 

Pu 



5r xd x - 



V 



S Bso yd 3 x + O{r- 4 )(A0) 



z 

1 

r 

I- 1 - 

r 

1 



= 1- 



Js By xd 3 x-^ Js By yd 3 x + 0{r 
J S B ^zd 3 x + 0{r~ 4 ), 
S x x i d 3 x + 0(r- 3 ), 
S^d 3 x + 0(r~ 3 ), 
S a ^d 3 x + 0{r- 3 ), 



x 
x 



— r / Sp x xd 3 x — — 



J S Px yd 3 x + 0{i 



1 )(41) 
(42) 
(43) 
(44) 
(45) 

),(46) 



P Z = T 



Sp xd 3 x — - 

v r 

S Pz zd 3 x + 0(r~ 4 ), 
S v x l d 3 x + 0{r~ 3 ). 



S Py yd 3 x + 0(r- 4 ),(A7) 



(48) 
(49) 



Here we briefly explain our order of computing 10 
Poisson-type equations. First, we solve Bi and x since 
the source terms do not depend on other unknown space- 
time variables. Once we derive the symmetric tensor A^ 
from Bi and x, we solve ip iteratively till the convergence. 
After that we solve aip iteratively till the convergence. 
Finally, we solve Pi and x- We use Modified Incomplete 
Cholesk y decomposition Conjugate Gra dient (MICCG) 
method ijMurata. Natori fc Karakll990j) to solve elliptic 
equations. 

Finally, we should maintain our grid resolution at the 
central core of the collapsing star. The nested grid is one 
of the most appropriate way to handle this problem (see, 
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e.g.. lRufferJl 9921. Our method is a mimic version of the 
nested grid to handle the collapsing s tar problem, which 
is bas ed on the rigridding method of iShibata fc Shapiro! 
( 2002). For the post homologous collapse of the SMS, 
there are two different timcscales in the collapsing star; 
one is the collapsing timescale at the center and the other 
is the one at the envelope, which i s sufficiently longer 
than t hat at the center. Therefore, IShibata fc Shapirol 
( 2002) add the grid number at the intermediate stage of 
the collapse in order to control both regions, center and 
envelope of the star. Since we perform the collapse of 
an SMS in 3D, we concentrate on the central core of the 
collapse due to the limitation of computational resource. 
In order to maintain spatial grid resolution especially at 
the center, we shrink the total grid to one half as the core 
radius of the star becomes half from the beginning of the 
present grid resolution. Although this method changes 
the boundary condition when we "zoom in" the compu- 
tational domain, it is approximately good when the dif- 
ference of the gravitational mass between the end of the 
previous computational grid and the beginning is small. 
In fact, the difference of the gravitational mass between 
every two computational domains is less than w 10~ 4 , 
and therefore this "zoom-in" method should behave as a 
reasonable approximation for a collapse of differentially 
rotating SMSs. 

We monitor the gravitational mass M and the angular 
momentum J 

M=-^-<f \7 l ipdSi 



= / [[(p + pe 



+ P){au t y - P] if; 



1 

16tt 



d 3 x, 



~j [xKl-yKi^dSi 



(xJ y - yJ x )ip & d 3 x 



(50) 



(51) 



during the evolution. In all cases reported in § 01 the 
gravitational mass was conserved up to ~ 0.1%, and the 
angular momentum up to ~ 1% of their initial values. 

We also compute proper mass M p , rotational kinetic 
energy T, gravitational binding energy W in equilibrium 
as 

Af„: 



T = 




e)d 3 x, (52) 
(53) 

W^M P +T - M, (54) 
where fi the angular velocity of the star. 

Since we use a polytropic equation of state at t = 0, 
it is convenient to rescale all quantiti es with respect to 
k JCook. Shapiro fc Teiikolskvl Il99l . Since K n / 2 has 
dimensions of length, we introduce the following nondi- 
mensional variables 

t = n- n / 2 t, x = n- n / 2 x, 



y = K n ? 2 y, 



z = k- ,1 / 2 z, Q = K n / 2 Q, M = K- n ' 2 M, (55) 




Fig. 1. — Comparison between numerical and analytical results 
in a one-dimensional relativistic wall shock problem at t/M = 1.0 
(where the fluid flow is aligned with the x-axis). Solid and dashed 
lines represent analytic and numerical results, respectively. For 
this simulation we chose T = 1.33, = 1.00 X 10~ 2 with a grid 
space Sx = 5 X 10 -3 and v = 0.200. 



Henceforth, we adopt nondimensional quantities, but 
omit the bars for convenience (equivalently, we set n = 
!)■ 

3. CODE TESTS 

First, we demonstrate ID relativistic wall shock prob- 
lem to check whether our code has an ability to treat 
shock. We set up the form of art ificial pressure viscosity 
ijHawlev. Smarr. fc WilsonlllQSl as 



Cvis p* 

o, 



Te)(Sv) 2 , for Sv < 0; 

for Sv > 0, 



(56) 



R = k- ,1 / 2 R, J = K~ n J. 



where Sv = 2SxdiV z , Sx(= Ax = Ay = Az) is the 
local grid spacing and where C V i S is the dimensionless 
parameter. When evolving the matter equations we 
limit the stepsize At by the following Courant condi- 
tion (At = min(0.3, Cd yn / V Pmax) Ax) , where the second 
term represents dynamical time of a collapsing star. We 
choose the dimensionless parameter Cdyn ~ 0.01. 

We have tested the ability of our code to resolve shocks 
by performing a wall-shock problem, in which two phases 
of a fluid collide. In Fig. ^we compare numerical results 
with the analytic solution for initial velocities that are 
similar to those found in our simulations below. With 
C v is = 3 we find good agreement, and set this value to 
simulate the catastrophic collapse in <2J 

Next, we demonstrate a spherical dust collapse in Fig|3 
Note that our conformally flat approximation retains all 
the nonlinear terms to maintain the exact dynamics for 
a spherical spacetime. We compare ID and 3D results 
to check whether our 3D code has an ability to repro- 
duce a spherical dust collapse. Note that we construct 
our ID code following Wilson's approach l)Wilsonlll979() 
described in Appendix El We choose the grid size of ID 
computation as 5,000, while of 3D one as (101 x 101x51). 
We find a very good agreement of the central lapse be- 
tween our ID result and our 3D conformally flat one 
within the error of 1%. We terminate the integration 
of our 3D code at t/M — 18.6, since the convergence of 
the iteration process in our elliptic solver becomes signif- 
icantly slow. 

Then, we demonstrate the collapse of a spherical star 
in Fig. |3 (see Tabled for an equilibrium profile). Since a 
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Fig. 2. — Comparison of the central lapse between three- 
dimensional and one-dimensional results in Oppenhymer-Snyder 
collapse. We start the collapse of a spherical dust from R = 
4M. Solid and dashed lines represent the central lapse of three- 
dimensional and one-dimensional results, respectively. 



TABLE 1 

Parameters for the initial spherical 
equilibrium sms 



Po 




M c 


Rc/M 


3.80 x 10~ 6 


2.49 x 10 2 


3.83 


65.0 



a Maximum rest-mass density 

b Equatorial circumferential radius 

c Gravitational mass 



radially unstable spherical star only collapses to form a 
spherical BH promptly, we should follow the collapse of a 
spherical star to check whether our code has an ability to 
follow the collapse and to find a signal of BH formation. 
We set up the same grid size as we did for a spherical dust 
collapse, covering the star with 81 grid points across the 
diameter. The central lapse of the star decreases mono- 
tonically, which means a collapse of a spherical star di- 
rectly forms a spherical BH. Although we terminate our 
integration around a c ~ 0.1, the fact that ID compu- 
tation of the spherical star collapse finds an apparent 
horizon at a c ~ 0.11 indicates that a BH might form. 

Finally, we check the stability of a uniformly rotating 
star whether our code has an ability to determine the 
radial stability of a star. Since we determine the radial 
stability of a differentially rotating star in 21 we should 
check the sensitivity of our code to determine the criti- 
cal onset of radial instability of uniformly rotating stars 
by comparing the results derived from the turning point 
method. 

To assess the ability of our code to distinguish sta- 
ble stars from unstable ones with rotation, we consider 
an equilibrium sequence of uniformly rotating stars of 
the fixed angular momentum J (J/M 2 = 0.644 at the 
turning point). While the turning point criterion strictly 
identifies the onset of secular instability the point of on- 
set of dynamical instability nearly coincides with the sec- 
ular instability point. We adopt the polytropic index of 
n = 2.96, which is regarded as radiation pressure domi- 
nant, SMS sequence, and grid resolution as in the spher- 
ical simulations reported above. We decrease the initial 
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Fig. 3. — Central lapse of the spherical collapse (see Table. 0. 
Monotonic decrease of central lapse indicates a prompt collapse to 
a spherical BH. 
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Fig. 4. — Probing the dynamical stability of a rotating SMS with 
n = 2.96, J = 10. Here, p c is the central density of the equilibrium 
rotating star. Filled circles and crosses represent unstable stars, 
while open circles represent stable stars according to our dynam- 
ical calculation. A cross indicates that the star is actually stable 
analytically according to the turning point criterion. The radii of 
the 5 marked stars are R/M = 254, 421, 539, 708, and 1579, where 
the sequence starts at the right side of the figure at the highest 
central densities. Note that the solid line shows a constant J se- 
quence with J = 10, while the dashed line represents the spherical 
equilibrium sequence. With the adopted grid resolution, our code 
can distinguish stable stars from unstable ones within 0.3% of the 
maximum gravitational mass. 



pressure to induce the collapse (k — ► 0.99k). 

Figure 0] summarizes our dynamical stability analysis 
for the rotating SMS. We conclude that with the adopted 
grid resolution, our code can distinguish stable rotating 
stars from unstable ones within 0.3% of the maximum 
gravitational mass. Figure shows the evolution of the 
central density for stable and unstable rotating stars. 

4. COLLAPSE OF A DIFFERENTIALLY ROTATING 
SUPERMASSIVE STAR 

Here we explain our initial data sets for collapsing 
stars. Since we are interested in a collapsing star of 
J/M 2 ~ 1, we have three requirements to construct dif- 
ferentially rotating equilibrium stars. 
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Fig. 5. — Evolution of the central densities of the stars plotted 
in Fig. [3] Curves are drawn for stars which are unstable both 
numerically and according to the turning point criterion (solid), 
unstable numerically but stable according to the turning point cri- 
terion (dashed), and stable both numerically and according to the 
turning point criterion (dash-dotted). 



The first requirement is to construct radially un- 
stable differentially rotating stars. The critical T to 
the onset of radial instability for slow rotation and 
weak gravitational field is described analytically as 
llChandrasekhar fc Lebovit^ 119681 iShapiro fc Teukolskvl 
1985) 



_ 4 M 2 n 2 i 
T crit _ - + 2.25— -- — , 



(57) 



where / is the inertia momenta of the star. Note that 
relativistic gravitation unstabilizes the star, while rota- 
tion, which produces centrifugal force, stabilizes the star. 
Although we take into account the differential rotation in 
the equilibrium star, this criterion (eq. |57|) may indicate 
the appropriate direction to choose a parameter sets of 
the initial condition of the rotating collapsing star. From 
eq. H57|) , we should at least choose soft equation of state 
to induce the collapse and set aiin = 3 polytropic star, 
an SMS sequence of the star. 

The second requirement is to construct a star which 
holds J/M 2 1. From the critical onset of radial insta- 
bility in the equilibrium star, a uniformly rotating SMS 
takes the maximum of J/M 2 ~ 0.9 when R/M ~ 600 
l|Baumgarte fc Shapiro! Il999t lShibatall2004|) . To verify 
the nature of cosmic censor, we should at least go be- 
yond J/M 2 > 1. The main restriction to hold large 
J/M 2 comes from the mass shedding limit of the star. 
Therefore, to construct a star of J/M 2 > 1 in radially 
unstable branch, differential rotation of the star is re- 
quired in SMS sequence. 

The last requirement is high degree of differential ro- 
tation. It is indeed one path to construct high degree of 
differential rotation in quasi-static evolution, and reach 
the onset of radial instability. 

To summarize, we need soft equation of state and high 
degree of differential rotation, and choose ann — 3 poly- 
tropic star (SMS sequences) and the degree of differential 
rotation as £l c /tt cq « 10 where we define the rotation 
profile as 

u t U{p = A 2 {Vl c -n). (58) 



In the Newtonian limit (it* 
tion law can be written as 



1, u u 



A 2 n r 



A 2 '' 



zu 2 il), this rota- 



(59) 



where A is the degree of differential rotation, vj is the 
cylindrical radius of the star. Since A has a dimension 
of length, we normalize it with a proper equatorial ra- 
dius R e , (A = R e A). Hereafter we choose A = 1/3 to 
construct relatively a high degree of differential rotation. 
We briefly summarize our method to construct relativis- 
tic rotating equilibrium stars in Appendix iBl 

We summarize the parameters of our differentially ro- 
tating stars at initial in Table |2 We slightly perturb our 
initial equilibrium state according to 



p = ^(equilibrium) ( ]_ + #(1) X + V + fl(2) ; 

\ R c 



y 



R 2 



,(60) 



where £W = §( 2 ) = 10~ 3 . We install m = 1 and m = 2 
density perturbation to provide the seed for one-armed 
spiral and bar formation, if the physical situation should 
lead to unstable growth. We adopt a grid size (201 x 
201 x 61), so that the star is initially covered by 161 
points across the equatorial diameter. We evolve the 
rotating SMS up to the point at which the conformally 
flat approximation breaks down. 

Figure [S] shows the results of radial stability in 4 dif- 
ferentially rotating stars. Since we do not have a tool 
to determine radial stability in differentially rotating 
stars from their equilibrium states (but for determin- 
ing th£_crit£non_of_s^c^^ stars, 
see lBonazzola. Frieben fc Gourgoulhonl l)1998j) ). the evo- 
lution is necessary to determine its stability. The crite- 
rion to determine the radial stability is as follows. When 
the central density of the star grows exponentially within 
a few dynamical time, we determine the star radially un- 
stable. On the other hand, when the central density of 
the star oscillates around its equilibrium state, we deter- 
mine the star radially stable. From this criterion, Model I 
and II have an exponential growth of the central density 
that means radially unstable, while Model III and IV 
have maximum at several dynamical times that means 
radially stable. 

We also show the evolution of central lapse in Fig. [7| 
The rapid decrease of a c below 0.3 indicates that a BH 
is likely to form. Model II shows that the central lapse 
monotonically decreases from ~ 0.9 to ~ 0.3. This figure 
shows that we can follow the collapse from the regime 
of Newtonian gravity (a c ~ 0.9) to that of relativistic 
gravity (a c ~ 0.3). On the other hand Model III shows 
that the central lapse oscillates around its equilibrium 
state, and that also means the star is radially stable. 

We show the final density snapshots of the stars in the 
equatorial plane (Fig. |SJ and in the meridional plane 
(Fig. Even in the final snapshots of radially unstable 
stars, the collapse is almost axisymmetric. Also, from the 
snapshots in the meridional plane, the material around 
the rotational axes collapses faster than the surrounding 
material due to the strong nonlinear gravitational field. 

Let us now focus on the collapsing star (Model II) to 
probe the final outcome. We locally define the cylindri- 
cal rest mass m, angular momentum j, specific angular 
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TABLE 2 

Parameters for the initial differentially rotating equilibrium 

SMSs 



Parameter 


Model I 


Model II 


Model III 


Model IV 


Rp/Re a 


0.600 


0.575 


0.550 


0.500 


.max 

Co 


3.38 x 10~ 6 


3.38 x lO" 6 


3.38 x lO" 6 


3.38 x lO" 6 


Rc 


3.06 x 10 2 


3.12 x 10 2 


3.18 x 10 2 


3.35 x 10 2 


n c b 


1.45 x 10~ 3 


1.49 x 10~ 3 


1.53 x 10~ 3 


1.59 x 10~ 3 




1.38 x 10~ 4 


1.42 x 10~ 4 


1.45 x 10~ 4 


1.51 x 10~ 4 


M 


4.69 


4.78 


4.88 


5.10 


J d 


2.13 x 10 1 


2.31 x 10 1 


2.51 x 10 1 


2.96 x 10 1 


T/W e 


6.47 x 10~ 2 


7.02 x 10~ 2 


7.60 x 10~ 2 


8.80 x 10~ 2 


J/M 2 


0.97 


1.01 


1.05 


1.14 


Rc/M 


6.53 x 10 1 


6.52 x 10 1 


6.52 x 10 1 


6.57 x 10 1 


Stability 


Unstable 


Unstable 


Stable 


Stable 



a Ratio of the polar proper radius to the equatorial proper radius 

^Maximum rest-mass density 

c Equatorial angular velocity at the surface 

^Angular momentum 

°Ratio of the rotational kinetic energy to the gravitational binding energy 




Fig. 6. — Evolution of the central density of 4 differentially 
rotating stars. Solid, dotted, dashed and dash-dotted line denotes 
Model I, II, III and IV (see Table 




Fig. 7. — Evolution of the central lapse of 2 differentially rotating 
stars. Solid and dashed line denotes Model II and III, respectively 
(see Table 0. 



momentum j s , Keplerian angular velocity fix as 

p CO rVJ 

m — 4:TT dz dwwp*, 
Jo Jo 



j = 4?r 



dz 



dwwhu 




j s =4:TT / dz dwwp*hu v , 



(61) 
(62) 
(63) 
(64) 



to indicate the transport of angular momentum and the 
distribution of j/m 2 . We assume the axisymmetric col- 
lapse to investigate the final outcome, since the merid- 
ional density snapshots (Fig. [§)) behave almost axisym- 
metry. 

Figurc lTOl shows the concentration of mass density pro- 
file during the collapse. As the collapse goes on, the in- 
crease ratio of the cylindrical mass at the core of the star 
becomes high. In fact, 75% of the cylindrical mass is in- 
side the radius of r < 2M at t = 5.64td yn - We also show 

the mean radius, defined as r m = (J dvp*-uc 2 )/M^, dur- 
ing the collapse in Fig. ^] Note that M* is the rest mass. 
We find that the mean radius monotonically decreases, 
which means that the collapse of the central core is co- 
herent. 

Figure 1121 shows the variation of angular velocity pro- 
file during the collapse. Since the cylindrical rest mass 
density shows a coherent collapse, the degree of differ- 
ential rotation significantly increases at the core during 
the collapse and approaches to the "Keplerian" angular 
velocity defined in eq. I|64|l . 

We show the specific angular momentum distribution 
during the collapse in Fig. In a global sense, J/M 2 
is conserved so that a final BH should violate cosmic 
censor, if it forms. However, we find that most of the 
matter collapses to form a BH, while not a few amount 
of angular momentum stays at the surface area of the 
star that prevents forming a BH of J/M 2 > 1. 

We show the distribution of j/m 2 in Fig. ^] during 
the collapse to verify cosmic censor. We find that mass 
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Fig. 8. — Final density contour in the equatorial plane of 4 differentially rotating stars. Model I, II, III, IV is plotted at the parameter 
(t/t dyu , p^ ax ) = (4.38, 1.37 X 1CT 3 ), (5.62, 2.16 X 10" 2 ), (5.61, 9.13 X 1CT 6 ), (5.60, 5.33 X 10" 6 ), respectively. The contour lines denote 
densities p* = p^ ax X lo-° 267 ( 16 - i ) (j = 1, . . ■ , 15). 



density collapses first in the central part and angular mo- 
mentum remains in the surface of the star that prevents 
forming a BH of J/M 2 > 1. Note that the distribution of 
j /to 2 is only an indicator to interpret the physical cause 
to prevent BH formation. Although we should define the 
total gravitational mass locally to discuss the local dis- 
tribution of J/M 2 , there is no knowing how to define it. 
Therefore we use the local rest mass instead. 

5. DISCUSSION 

We investigate the collapse of differentially rotating 
SMSs by means of hydrodynamic simulations in confor- 
mally flat approximation in general relativity. We start 
our collapse around the onset of radially instability at 
R/M ~ 65 to the point where conformally flat approxi- 
mation breaks down. 

We find that cosmic censor even holds for a gravita- 
tional collapse of a radially unstable differentially rotat- 
ing equilibrium SMS of J/M 2 > 1. The main reason to 



prevent formation of a BH of J/M 2 > 1 is that quite a 
large amount of angular momentum stays at the surface, 
not core bounce nor shock propagation, bar formation in 
our model. Note that even a thin disk near the surface 
of the star can hold relatively a large amount of angular 
momentum if the radius is large. The abo ve conclusion is 
supported by Sekiguchi & Shibata ( 2004) in full general 
relativistic simulation that the criterion of BH formation 
is determined at the central q c (= j /m 2 \ tu ^o)- In fact, 
the central q c of Model H is q c ss 0.89, which satisfies 
their criterion. 

The collapse of a differentially rotating, radially un- 
stable SMS is coherent and likely leads to formation 
of an SMBH. This situation is quite s imilar to the col- 
lapse of a uniformly rota ting SMS (e.g.. lSaiio et alJl2002t 
iShibata fc Shapiroll2002]) . since T/W « 0.08 seems still 
not sufficient to change the path of the collapse from 
the one of uniform rotation. However, this final out- 
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come may depend on the equation of state of the star. 
lLoeb fc Rasiol ()1994fl treated the isothermal (T = 1) 
collapse of initially homogeneous, uniformly rotating, 
low entropy clouds via smooth particle hydrodynamics 
(SPH) simulations. They found considerable fragmen- 
tation into dense clumps, and disk formation containing 
~ 5% of the mass. They concluded that a seed BH will 
form at the center and that it li kely will grow g radually 
by accretion in our model. Also Shibata (2003|) treated 
the collapse of uniformly rotating polytropic star from 
the critical onset in the range of T w 1.5 — 2.5 and found 
that the mass of the disk is less than 10~ 3 of the gravi- 
tational mass. 

We cannot find any evidence of bar formation nor sig- 
nificant disk formation from the rotating collapse prior 
to BH formation. The phenomenon of no bar formation 
also comes from the fact that mass density collapses first 



to form a BH. In such case, T/W cannot scale in due 
to the growth of the degree of differential rotation, and 
as a fact the star of T/W cannot reach the dynamical 
instability point of ~ 0.27. Since the m = 1 dynamical 
instability takes place when the st ar has a toroidal struc- 
ture and soft equation of state llCentrella et alJ 120011 
iSaiio. Baumgarte. fc Shaniroll2003j) . we cannot find any 
evidence of toroidal structure at the time we stop our 
integration. 

The rotating SMS collapse is a promising source of 
burst gravitational waves and of quasi- normal mode ring- 
ing waves. The characteristic frequency for burst (/burst)) 
quasi- normal mode ringing (/qnm), and the wave am- 
plitude for burst (ft-burst), quasi- normal mode ringing 
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Fig. 10. — Mass density profile in the x-axes for Model II. Note 
that M t is the rest mass. Solid and dashed line denotes the profile 
at t = and t = 5.64 t Ayn , respectively. 
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Fig. 13. — Specific angular momentum profile as a function of 
cylindrical mass for Model II. Solid and dashed line denotes the 
profile at t = and t = 5.64 t dyn , respectively. 





FIG. 11.- Mean radius of the star during the collapse for Model M ° de L! L S ° Ud and daS , hed line den ° teS th ° pr ° filc at * = and 

r t = 5.64 tdym respectively. 




Fig. 12. — Angular velocity profile in the x-axes for Model 
II. Solid and dashed lines denote ang ular velocity of the star and 
Keplerian angular velocity (see eq. 1641 for its definition). The 
top panel shows the snapshots of t = while the bottom panel of 
t ~ 5.64 t dyn . 
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(66) 
(67) 

1/2 

(68) 



where d is the distance from the observer and AEqw 
is the total radiated energy. We set R/M = 1, a char- 
acteristic mean radius during BH formation. Since the 
main targets of LISA are gravitational radiation sources 
between 10~ 4 and 10 _1 Hz, it is possible that LISA can 
search for the burst and quasi-normal ringing waves ac- 
companying rotating SMS collapse and formation of an 
SMBH. 



(/iqnm)) are ijSaiio et al.ll2002|) 

/b „„.3x 1 0-(^fe)(f) 3,2 |Hzl , (65) 
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APPENDIX 

RELATIVISTIC HYDRODYNAMICS IN SPHERICALLY SYMMETRIC SPACETIME 

Here we briefly explain relativistic hydrod ynamics in spherically symmetric spacetime (Wilson [19791 
IShapiro fc TeukolskvllT980l IBowers k WilsoiJl99l|) that we use in § |21 We can describe the line element of spherically 
symmetric spacetime in the isotropic coordinate as 

ds 2 = (-a 2 + ^(3 2 )dt 2 + 2ip A (3dtdr + ^{dr 2 + r 2 dil), (Al) 

where a is the lapse, f3 is the shift, if) is the conformal factor. 

We use the maximal slicing condition K — 0, which derives the following relation to the components of an extrinsic 
curvature as 

K = -2K e e {= -2K*). (A2) 

As a consequence, we only need to consider K^, to construct an extrinsic curvature. 
The momentum constraint derives the equation of the extrinsic curvature as 

^(riJj 2 K r r ) = (rTp 2 ) 3 8nJ r . (A3) 

The integration form of eq. I|A3J) is 

K = 7^)3 fir^Jrdr. (A4) 
The restriction of the spatial metric to be conformally flat requires the following equation 

2aK r r = \rd r l-\ , (A5) 



3 

which gives us an appropriate boundary condition for the extrinsic curvature as 0(r -3 ). Also the integration form of 



eq. 



P = | r / ^Ldr. (A6) 
2 ./„ r 



The Hamiltonian constraint and the trace of the evolution equation guide 

^d r {r 2 8A) = -2WVh - ^(K) 2 , (A7) 
r A 16 

\d r \r 2 d r {a^)\ = 2W 5 (p H + 2S) + ^ a ^{K r r ) 2 , (A8) 
r 16 

with the boundary conditions 

V^ = l + ^- + (r- 3 ), (A9) 
2r 

M 

aip = l + — +o(r" 3 ), (A10) 
r 

at the grid edge and 

d r a = d r i) = 0, (All) 

at the center (r = 0). 

Therefore, we can determine 4 unknown variables (K£, ip, at/), (3) with 4 equations (eqs. |A4) . |A6j . |A7j . |A8j ). 3 
Matter equations can be written as 

d(r 2 p*) d{r 2 p*v r ) 



dt dr 
9(r 2 e*) d(r 2 e«v r ) 
dt dr 



: 0, (A12) 

0, (A13) 



d{r 2 p*u r ) d{r 2 p*u r v r ) 2 6 _ t 

— 1 = -r QfP r - p*au a, r 

ot or 

+r 2 Pr u r (3\ + r 2 ^4 L ^ r . (A14) 



3 We choose A£ = ip®K£ instead of K£ from a computational point of view (See i|2.ijl '). 
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We have not taken into account artificial viscosity, since shock does not seem to play an important role in a spherical 
dust collapse and the collapse of a spherical star. In §0 we adopt our ID relativistic hydrodynamic code to a spherical 
dust collapse by setting e* and <C equilibrium variables of them in typical polytropic index. 

We briefly mention a method to compute an apparent horizon in ID (Sha piro fe Teukolskvl Il980t 
iNakamura. Oohara. fc Koiimal ll987). An apparent horizon is defined as an outermost margin ally trapped surface 
whose future-directed outgoing null geodesies have zero expansion O (Hawkin gs fc Ellis! 119731) . The outgoing null 
vector can be expressed as k^ — (n^ 1 + s' 1 )/-^, where n M is unit normal of the spacelike hypersurface defined as 
n 1 * =(—a, 0, 0, 0), s M is the out-directed spacelike vector orthogonal to the surface. The projection tensor /i M „ and 



the extrinsic curvature Kij are described as 



The expansion is then written as 



e 



K. 



(A15) 
(A16) 

(A17) 



For spherically symmetric spacetime, we can choose the spacelike vector s l as s 1 = (ip 2 ,0,0), and then the radius 
tah of an apparent horizon satisfies the following equation, 



1 



2r % + \t^K\ 
ip 2 



0. 



(A18) 



Note that an event horizon always lies outside an apparent horizon, and both horizons coincide with each other 
when the spacetime is s tationary. Also an event horizon always forms before the time that an apparent horizon does 
ijHawkings fc Ellisl[T97^l . 



CONSTRUCTION OF RELATIVISTIC ROTATING EQUILIBRIUM STARS 

Here we summarize our method to construct rotating rela t ivistic equ ilibrium stars, which is based on 
iKomatsu. Eriguchi. fc Hachisul lfl989j) : iCook. Shapiro, fc Teukolskvl lfl"992t [l994) (see also lStergiouias1l2003l for a his- 
torical review). 

First we focus on how to solve the relativistic Euler equation. The equation can be described in axisymmetric 
spacetime as 



-T f + zfuMj = 0, 

h ir 



where h = (1 + e + P/po) is a specific enthalpy. We also assume specific type of rotation law as 



U U u 



(Bl) 



(B2) 



where VL C is the central angular velocity, in order to integrate eq. IjBljl . The Bernoulli's equation is driven by integrating 
the relativistic Euler equation (|B1|I as 

H - K + R= C, (B3) 



where 



H = 
K = 
R = 



h 

du 



ln[l + (n + l)q], 



1 



—r- = \nu f = — \u[a 2 - ip 4 [(/3 x - n y ) 2 + (f3 v + n x )% 



-A 2 (fl c -fl) 2 



(B4) 
(B5) 
(B6) 



where q = P/po 4 - Note that we have already assumed polytropic equation of state P = Pq wher e r = 1 + l/n, n is a 
polytropic index. Therefore, we can describe the matter distribution equation by using eq. I|B3() and the rotation law 
(eq. [H3) as 



1 



n + 1 
1 



Cexp[A 2 (n c ~n) 2 /2] 



n 



1 



Cexp[A 2 (n c 



- 1 

n) 2 /2] 



A 2 (fl c - f2) =u t u v 



i> 4 [x(l3 v + fix) - y{(3 x - fly)} 
a 2 - iP 4 [{f3 x - fly) 2 + (0y + fix) 2 } ■ 



- 1 



(B7) 
(B8) 



We only use q = P/po in this section. 
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Note that we assume conformally flat spacetime to derive the final part of the equality in eqs. i|B7|) and JB8I) . 

From the computational point of view, we introduce nondimensional quantities rescaled in the equatorial proper 
radius of the star (i? e ) as 

x = x/R e , y = y/R e , z = z/R e , 

{l = R e n, A = A/R e . A = R 2 e A. [ ' 

We also rescale the lapse and conformal factor using the nature of scale free in Newtonian gravity as 

a = a 1 /^> = V^ 1/(2/i? ' ) - (BIO) 

To determine the matter distribution and rotation profile of the next iteration step, we have to determine five 
unknown variables (R e , C, fl c , f2 e , fi max ) from five equations. Therefore, we evaluate eqs. (|BT|) and (|B8|) at three 
locations, the point of the maximum rest mass density, that of polar surface, and that of the equatorial surface. Note 
that eq. i|B8(l becomes a n identity at pola r surface of the star, and thus we have five equations to be solved. We use 
New ton- Rapson method l)Press e t al. 1992]) to solve these equations. Once we determine 5 unknown variables, we solve 
cq. (|B8|) using Newton-Rapson method to determine the rotation profile. After that we can determine the matter 
distribution (eq. |B7j h 

Next we briefly summarize the gravitational field equations normalized by proper radius as follows. 

AB, = SirRl^J, = AttSb, , (Bll) 
Ax = -SvRl^JrX 1 = 4tt5 x , (B12) 

Aip = ~2ttR1^ 5 Pr - \^- 7 A t] A^ = 4^, (B13) 
8 

A(m/0 =2-KRlail){p a + 25) + 7 -a^ 7 A %3 A 11 = 4tt5 q ^, (B14) 

8 

APi = AnRlaJi = 4nS Pi , (B15) 

Ar) = -kirRlaJiX i = 4irS n . (B16) 

We choose the same boundary conditions as we choose in our evolution code (see § EJ ■ Note that we construct the 
star in 3D using octant symmetry to adopt our initial data to our 3D evolu tio n cod e smoothly. 

To summarize, we first solve the gravitational field equations (eqs . IbTT] - |B16p . Next, we determine 5 unknown 
quantities (i? e , C, fl c , O e , £! max ) from 5 equations (eqs. |B7) and |B8j ) using Newton-Rapson method. Finally we 
determine the matter distribution and rotation profile for the next iteration step. We continue this iteration cycle till 
R e converges; i.e. R e goes bellow the error of 10~ 5 . We also check the error rate of physical quantities such as M and 
J and find that they are in fact below the error rate of 10~ 4 . 
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